A Real-Time Evaluation Algorithm for Noncontact Heart Rate Variability Monitoring

Noncontact vital sign monitoring based on radar has attracted great interest in many fields. Heart Rate Variability (HRV), which measures the fluctuation of heartbeat intervals, has been considered as an important indicator for general health evaluation. This paper proposes a new algorithm for HRV monitoring in which frequency-modulated continuous-wave (FMCW) radar is used to separate echo signals from different distances, and the beamforming technique is adopted to improve signal quality. After the phase reflecting the chest wall motion is demodulated, the acceleration is calculated to enhance the heartbeat and suppress the impact of respiration. The time interval of each heartbeat is estimated based on the smoothed acceleration waveform. Finally, a joint optimization algorithm was developed and is used to precisely segment the acceleration signal for analyzing HRV. Experimental results from 10 participants show the potential of the proposed algorithm for obtaining a noncontact HRV estimation with high accuracy. The proposed algorithm can measure the interbeat interval (IBI) with a root mean square error (RMSE) of 14.9 ms and accurately estimate HRV parameters with an RMSE of 3.24 ms for MEAN (the average value of the IBI), 4.91 ms for the standard deviation of normal to normal (SDNN), and 9.10 ms for the root mean square of successive differences (RMSSD). These results demonstrate the effectiveness and feasibility of the proposed method in emotion recognition, sleep monitoring, and heart disease diagnosis.


Introduction
Heart Rate Variability (HRV) describes the variation in the time interval between successive heartbeats. It is generally considered to be the result of an interaction between the heart and the brain, which is also called neuro-cardiac function. HRV is controlled by the autonomic nervous system (ANS), including the parasympathetic nervous system (PNS) and sympathetic nervous system (SNS). In recent years, research on HRV has attracted lots of interest from the fields of instrumentation [1], signal processing [2], and healthcare [3]. Studies have confirmed that HRV indicators can be applied to sleep monitoring [4] and the diagnosis of some cardiovascular diseases [5]. They can also be used to evaluate people's mental stress [6] and can even further help home buyers choose residences [7].

Background
Electrocardiography (ECG) and photoplethysmography (PPG), which reflect the change in body surface electrical potentials and blood volume fluctuations in a super-ficial body location, respectively, are still the mainstream monitoring methods for HRV [8]. However, as contact measurement methods, they both have many limitations in practice. Their lead wires will hinder the user's physical movement. For electrocardiographs, many electrodes need to be attached to particular locations, which complicates the monitoring process, makes people feel uncomfortable, and puts higher requirements on operators. In addition, they do not apply to certain groups, such as infants, patients with skin burns, and individuals with sleep disorders [9][10][11][12][13].
In the past 50 years, contactless physiological monitoring has developed rapidly [14] and has largely overcome the disadvantages of contact equipment such as electrocardiographs. It can monitor respiration, heartbeat, or other physiological signals without directly contacting the body of the monitored subject and will not make them feel uncomfortable or interfere with their daily routine [15], which allows its use for special groups like infants. Monitoring based on millimeter-wave radar is considered a promising noncontact physiological monitoring method that can penetrate non-metallic obstacles, such as clothes or quilts on the surface of the monitored subject, to capture physiological signals of the human body [16]. It can be used to probe respiration disorders, such as obstructive sleep apnea (OSA) and sudden infant death syndrome (SIDS), as well as used in medical sleep labs and earthquake or fire search-and-rescue scenarios [17,18]. The radiofrequency (RF) signal is more robust to temperature changes or environmental thermal noise compared with infrared thermal imaging methods [12] and is better able to avoid insufficient image resolution, blind areas, or potential privacy problems compared with vision-based monitoring methods [19].
However, measuring the time interval between successive heartbeats and analyzing HRV remain significant challenges due to the weak amplitude of the heartbeat signal and interferences from respiration, trunk movement, and various noises. Most of the existing research in this field focuses on monitoring the respiration rate and heartbeat rate [20], while only a few studies focus on the extraction of HRV characteristics.

Related Works
In 2007, Massagram et al. successfully extracted the time interval between successive heartbeats from echo signals of continuous-wave (CW) radar, proving the feasibility of monitoring HRV with the RF method [21]. Early studies on this subject mainly chose CW radar as RF equipment [16,21]. In some of these studies, the peak corresponding to the heartbeat signal can be clearly distinguished in the demodulated chest movement signal. So, the heartbeat signals can be easily separated by using some filters or using classical signal processing analysis methods, such as autocorrelation, from which the time interval between successive heartbeats can be extracted so that HRV can be measured [22][23][24]. However, CW radar has a weak anti-interference ability. When other subjects appear in the monitoring range, the movement of one individual will affect the measurement and extraction of physiological signals from the other individual, leading to a large measurement error.
FMCW radar provides an important approach to solve this problem. The modulated frequency provides a unique range resolution that CW radar does not have, which enables the separation of the echo signals of objects at specific distances. As a result, FMCW has a strong anti-interference ability and the potential for monitoring multiple subjects simultaneously [25,26]. Effective algorithms for determining the distance between monitored subjects and radar were also developed [27][28][29]. When there are multiple individuals at the same distance from the radar, their vital signals can still be separated and monitored simultaneously with beamforming technologies [30][31][32].
In practice, the amplitude of the heartbeat waveform demodulated from the FMCW radar echo signal is weak: it is orders of magnitude smaller than the amplitude of the respiration waveform and almost buried in a combination of harmonic respiration signals and noise in the frequency spectrum of the phase sequence. Therefore, the design of the heartbeat signal extraction algorithm plays an important role in monitoring.
A bandpass filter is a common method to separate heartbeat signals from demodulated chest wall motion signals because it is simple to design and easy to implement [23,33,34]. However, it is possible to mistake the harmonics of respiration for heartbeat signals because of their similar amplitudes and frequencies. In addition, it is difficult to accurately extract the time interval of each heartbeat because it is difficult for the filter to eliminate the influence of respiratory harmonics. Most previous studies performed HRV analysis based on the average heart rate over a short time (e.g., 3 s) [33,34]. Empirical mode decomposition (EMD) is another commonly used heartbeat signal extraction algorithm [35][36][37] that decomposes the chest wall motion according to the time scale characteristics of the signal itself. However, this method has limitations, such as mode aliasing and end effects. In addition, selecting the components of the heartbeat signal from the intrinsic mode function (IMF) remains a technical challenge.
With the development of signal processing technology, some novel algorithms were proposed to extract heartbeat signals and achieved accurate results. Lv et al. used the envelope mid-line method to remove the respiration signal and proposed a stochastic resonance algorithm to enhance the amplitude of heartbeat signals [30]. The experimental results were highly consistent with those of PPG detection. The average accuracy of the heart rate value in all subjects reached 96.56%, and the relative error of SDNN was less than 6.53%. Xiong et al. [38] proposed a differential enhancement (DE) method, which uses differential operation to significantly enhance the heartbeat components, especially high-order heartbeat harmonics. Combined with the autocorrelation-based periodicity extraction technique, DE can locate the true heartbeat rate (HR). However, it is challenging to accurately extract the duration of each heartbeat because the waveform after the difference operation is complicated, and it is difficult to segment through peak detection. This problem is addressed by the method proposed by Zhao et al. [39] with the assumption that successive human heartbeats should have the same morphology. Hence, the corresponding heartbeat waveforms should have the same overall shapes, while they may stretch or compress due to different beat lengths. On this basis, the method transforms the above segmentation problem into an optimization problem and solves the optimal "template" while seeking the optimal segmentation. In their experimental results, the extracted time intervals of each heartbeat are within milliseconds of the ECG signal. However, the algorithm requires a lot of interpolation operations during the segmentation process, which makes it difficult to monitor the HRV index in real time.
In this paper, an HRV monitoring algorithm based on FMCW radar is proposed to address the issues of existing methods. It is based on a single-chip multiple-input multipleoutput (MIMO) radar sensor, zooms in on human reflections, and neglects reflections from other subjects by using the range resolution characteristics of FMCW radar and beamforming techniques. It was found that the change process of the respiration waveform is slow and steady, while a heartbeat involves rapid contraction of the muscles. So, the second-order difference is calculated to enhance the heartbeat and suppress the impact of respiration. An approximated estimate of the duration of each heartbeat can be obtained by calculating the short-time average power of the acceleration signal. Finally, a joint optimization algorithm was developed and is used to adjust the duration of each heartbeat to meet the accuracy requirements of HRV analysis. Compared with the algorithm in [39], the proposed method is ten times to dozens of times faster. Experiments are presented to verify the performance of the proposed algorithm in accurately extracting HRV characteristics and reducing computation. The proposed algorithm shows great potential for real-time evaluations.
This article is organized as follows: Section 1 introduces the background and the related works and discusses their advantages and limitations. Section 2 presents the proposed algorithm. Section 3 describes the experimental validation protocols, and the results are illustrated and discussed in Section 4. Finally, conclusions are drawn in Section 5. Figure 1 shows the main configuration of the proposed method, a real-time HRV index monitoring algorithm based on FMCW radar. The algorithm can effectively extract each heartbeat interval from the echo signal of FMCW radar and accurately estimate the HRV index.

FMCW Radar Signal Preprocessing
For FMCW radar, each transmitting antenna sends one linear frequency-modulated pulse (a chirp) at a time. Chest displacements due to breathing and heartbeat are measured through the phase modulation presented in the output of the mixer, delivering a signal that contains breathing and heartbeat effects. The mathematical expression of the intermediate frequency (IF) signal is denoted by where A is the amplitude of the IF signal, S is the slope of frequency over time, which is a constant, c is the speed of light, and λ is the wavelength of the impulse. R is the distance between the chest wall and the antenna. It is clear from (1) that the frequency f and phase ϕ of S IF (t) are proportional to the distance between the chest wall and antenna and can be shown as and ϕ = 4πR For each chirp, the IF signal is digitized by an analog-to-digital converter (ADC), producing N samples. The components of different frequencies in the output of the mixer can be separated by the fast Fourier transform (FFT). The component reflected by the human chest wall can be selected by (2) if the distance (donated by R 0 ) between the chest wall and antenna is known, from the phase of which the movement waveform of the chest wall can be obtained.
As shown in Figure 2a, the radar equipment in this research has three transmitting antennas and four receiving antennas. The distance between the adjacent transmitting and receiving antennas is λ and λ/2, respectively. The corresponding virtual antenna array is shown in Figure 2b. Based on the virtual antenna, the minimum variance distortionless response (MVDR) beamforming technology was adopted in our study to suppress interference reflected from other objects in the azimuth plane using the antennas numbered 1-8 in Figure 2b. Assuming that the azimuth relative to the antenna array of the monitored subject's heart is θ, the steering vector α(θ) can be written as where d = λ/2. The optimal weight vector can be formulated as where M is the number of virtual antennas, and R xx is the covariance matrix calculated by the receiving vector of the virtual antenna array.

Heartbeat Signal Enhancement
Although it is convenient to extract the heartbeat signal using a bandpass filter, some details of the heartbeat signal will also be filtered out because the heartbeat waveform is not strictly sinusoidal, and there will be residual respiratory harmonics, which is undesirable for performing HRV analysis.
We define the chest wall motion signal as R(n). Figure 3a shows an example of R(n), in which distinct breathing movements can be observed, but it is difficult to directly distinguish the vibrations caused by each heartbeat. Considering that the amplitude of vibration caused by the heartbeat in R(n) is small but intense, while the inhale-exhale motion is slow and smooth, the acceleration a(n) (i.e., the second derivative) of R(n) is calculated through (6) to suppress the respiration signal and enhance the heartbeat signal, where T s is the sampling interval. When the sampling frequency is 250 Hz, the corresponding T s is 4 ms.  Figure 3b shows the acceleration of the example shown in Figure 3a. The ECG signal sampled synchronously is shown in Figure 3c. It is clear that the heartbeat signal is enhanced and the breathing motion is dampened.
Note that while the periodicity of the heartbeat signal can be observed in the acceleration, it is still difficult to pick out each heartbeat cycle by using zero-crossing detection or peak detection because the waveform of the heartbeat in acceleration is not a simple sinusoidal signal, and the morphology of heartbeats in acceleration signals is unknown. For a single heartbeat in acceleration, the middle has a more violent fluctuation degree than the two ends. So, short-time average power P(n) is calculated to describe this feature and smooth the acceleration waveform. P(n) can be calculated as where L is half the length of the time window. If L is too small, the smoothing effect is not good, and if L is too large, the fluctuation of the waveform after smoothing is not obvious. Based on trail and error, the length of the time window in this paper is 0.4 s, which is about half of the heartbeat cycle. Figure 4 shows the acceleration signal before and after smoothing. The heartbeat signal is sinusoidal after smoothing, so it can be segmented easily. However, Figure 4 shows that many spikes in the smoothed waveform are not completely aligned with the ECG signal, which indicates that some details are lost when the acceleration signal is smoothed. Therefore, the smoothed waveform can only be used to roughly estimate the duration of each heartbeat but does not meet the accuracy requirements for HRV analysis.

Accurate Extraction of Heartbeat Intervals
The variation in the heartbeat time interval is usually a few milliseconds to tens of milliseconds, and the heartbeat interval extracted by the smoothed acceleration waveform does not meet this accuracy requirement. It is necessary to further improve the accuracy of heartbeat intervals for HRV analysis. Notice that although the waveforms of different single heartbeats are very complex, they are very similar in shape. Based on this, we can adjust the location of the cut-off point between two adjacent smoothed heartbeat waveforms. With these cut-off points, the template of a single heartbeat can be identified. This template can be used to further adjust the positions of cut-off points (denoted by S) between any two successive heartbeats. The solving process can be described as the following joint optimization problem: where S = s 1 s 2 s 3 · · · , s i is the cut-off point in the acceleration signal, a(s i−1 + 1 : s i ) is the sequence from the (s i−1 + 1)-th to the s i -th elements of acceleration, and ω(µ, p) indicates that the length of the template µ has been changed to p elements by cubic spline interpolation. The goal of the algorithm is to find the optimal segmentation S and template µ that minimize the sum of differences between each subsection and template, as shown in (8).
The optimization of the template and segmentation may require several iterations. In the iterative process, the segmentation and template are updated alternately, while the other one is fixed. During each iteration, the algorithm updates the template given the current segmentation through (9), and then it updates the segmentation with the updated template through (10). where The superscripts l and l + 1 indicate the number of iterations. For (9), it can be proved that the optimal template is given by (12), where m is the required length of the template, and n is the length of a(n). For (10), considering that the complexity increases exponentially with the length of a(n) when seeking the (l + 1)-th optimal segmentation S l+1 based on the l-th template µ l in the iteration, dynamic programming is adopted, as shown by (13), which reduces the complexity to a linear increase with the length of a(n).
where s i is a cut-off point, and s i−1 is a possible cut-off point, which is a specific range (e.g., 0.5-1.2 s) ahead of s i . However, if we select s i and s i−1 by using a traversal search, the amount of calculation will still be significant because a lot of interpolation operations are needed. In the proposed algorithm, it is assumed that the cut-off points in the current iteration can only appear next to the corresponding cut-off points that were obtained last time. Therefore, it is not necessary to include any two points within the range of the heartbeat cycle (e.g., 0.5-1.2 s) in the search for optimal segmentation. The real cut-off points can be found by just exploring the combination of points in the neighborhood of the estimated cut-off points from the smoothed a(n) or the last iteration, which can be written as where B is half the length of the neighborhood. The selection of the neighborhood size needs to take stability and computation into account. Too small a neighborhood size may screen the real cut-off point and lead to difficulties in converging. On the contrary, if the selected neighborhood is too large, the computation will increase when exploring the optimal segmentation, as the complexity will increase with the square of the neighborhood size. In this paper, B is selected as 20 ms. The value of B is 5 when the sampling frequency is 250 Hz. Figure 5 describes the algorithms for strengthening heartbeat signals and extracting heartbeat intervals. The algorithm can also be expressed by Algorithm 1. where

Algorithm Summary and Convergence Conditions
end for 16: end while 17: The time intervals between successive heartbeats are extracted from the last output S when the algorithm converges. It is necessary to define the convergence condition properly to decide whether to exit after each iteration is completed. It is noted that after several iterations, most of the cut-off points will become stable. For a certain cut-off point, its position will remain same. The convergence condition of the algorithm is defined as follows according to this characteristic.
where S l+1 j represents the j-th element of the set S l+1 , S l j is the corresponding element of the set S l , and #{S} represents the number of elements in S. If the difference between the two is less than a predefined threshold α 1 , the cut-off point is considered stable. When the proportion of stable cut-off points in all cut-off points reaches a threshold α 2 , the segmentation converges. In our study, considering that the time interval between two adjacent sampling points is 4 ms, the value of α 1 can be set to 5 ms to detect the stable points. The value of α 2 is set to 0.80.
The overall algorithm flowchart for HRV estimation based on radar echoes is shown in Figure 6. The chest wall motion signal R(n) is extracted through MVDR and range FFT based on R 0 and θ. The acceleration a(n) is calculated to enhance the heartbeat signal and smoothed by short-time average power. The positions of cut-off points are first estimated by the smoothed waveform and then adjusted by a joint optimization algorithm. Finally, the time intervals of successive heartbeats are obtained, and the HRV analysis is performed.

Experiments
A microwave radar evaluation module (EVM) produced by Texas Instruments (TI), IWR1843, was used to conduct the experiments. It has 3 transmitting antennas (Tx) and 4 receiving antennas (Rx). Each transmitting antenna can transmit 77-81 GHz RF signals. In the experiments, the start frequency of a chirp was set to 77 GHz, and the slope of frequency over time was 70 MHz/µs. Each chirp lasted 50 µs, during which 128 points were sampled. The fast time axis sampling frequency was set to 4 MHz, while the slow time axis sampling frequency was 100 Hz. Tx1, Tx2, and Tx3 transmitted a chirp as described above, in turn, during each frame. The patch antenna array had a peak gain of 10 dB, and the 3 dB beamwidth was approximately ±28 • . The mean power density was approximately 1 W/m 2 , which complies with the Exposure Limits in Council Recommendation 1999/519/EC. A DCA1000 EVM was used for radar data sampling, which is a capture card for interfacing with TI's 77 GHz IWR1843 EVM and enables users to stream the ADC data over Ethernet. The signal interface between the capture card and the IWR1843 EVM uses a 60-pin highdensity connector. Figure 7 shows the experimental setup for indoor lab testing and evaluation. The FMCW radar sensor was fixed 60 cm above the mattress by using a camera support. To measure the ground truth of IBIs, ECG monitoring was used. A total of 10 electrodes, including 4 attached to the extremities and 6 attached to the chest, corresponding to Wilson leads, were pasted on the body of the monitored subject. The experiment recruited 10 human subjects (8 male and 2 female; age: 28 ± 6 years; height: 175 ± 10 cm; weight: 58 ± 10 kg) with no respiratory or cardiovascular system diseases. The study was approved by the Medical Ethics Committee of Zhejiang Hospital (Approval Letter No.: AF/SC-06/04.2). The subjects were asked to lie down on the mattress and keep their torso still but breathe normally. The distance and azimuth of subjects' hearts were recorded at the beginning of the experiments, which were used to find the corresponding range bin and perform beamforming. When the subject was ready, the radar sensor was activated to start sampling first, followed by the ECG monitor. Although the sampled data of the two are not completely synchronized in time, data alignment can be easily achieved by observing the time intervals of each heartbeat and adjusting over a small range (e.g., 5 s). During the experiment, several sets of results with a length of 15 s were acquired for each subject. Four sets of data for each subject were retained after removing the poor-quality fragments. The data were exported to and processed by MATLAB. Note that the sampling frequency was changed to 250 Hz by interpolation before using the proposed algorithm to improve measurement accuracy.

Results
From 10 subjects, a total of 40 sets of data containing 694 heartbeats were obtained. Figure 8 shows all IBI measurement results. The scatter is plotted for each heartbeat, using the IBIs extracted from ECG as the horizontal axis and IBIs extracted from radar by the proposed algorithm as the vertical axis. The correlation coefficient between them is 0.9747. It can be seen that most of the points fall near the 45 • oblique line, indicating that the proposed algorithm has high measurement accuracy. The measurement error err of each IBI, defined in (16), was calculated to quantify the measurement accuracy. err(n) = IBI radar (n) − IBI ECG (n) (16) where IBI radar (n) and IBI ECG (n) are the interbeat interval of the n-th heartbeat extracted from radar and ECG, respectively.  Figure 9a describes the distribution of the measurement error. Most of the measurements have an error of ±2 ms, and the number of measurement points decreases for a large absolute value of error. The whole distribution is almost symmetric around zero error, indicating that the source of the measurement error is more likely a random error. Figure 9b shows the cumulative distribution function (CDF) of measurement results using the absolute value of the measurement error. Note that 36.46% of all the measurement points have an error within ±4 ms. When the allowed error range is enlarged to ±8 ms, over 60% (62.39%) of all points will meet the requirements. More than 80% of all measurement points are within 16 ms of the true value.
The root mean square error (RMSE) of IBIs is calculated by using (17). 17) where N is the total number of heartbeats recorded. The RMSE of the proposed algorithm is 14.9 ms for IBI measurements. The HRV features can be further obtained from IBIs. The three most widely used time-domain metrics are used to evaluate HRV. MEAN is the average value of IBIs, which can be calculated by using (18).
The standard deviation of normal to normal (SDNN) measures the standard deviation of all IBIs and can be calculated by using (19).
The root square of successive differences (RMSSD) measures the successive IBI changes and can be calculated as in (20).
The above three metrics were calculated for each set of sampling data, and the scatter of the results is drawn in the same form as Figure 8. Figure 10a shows the estimation results of MEAN. The correlation coefficient is 0.9976, indicating that the two have a strong correlation. The RMSE of MEAN in all 40 sets of data is 3.24 ms. Figure 10b,c show the estimation results of SDNN and RMSSD, respectively. Most of the points are above the oblique line, indicating that the corresponding HRV metrics calculated based on the RF signal are relatively large, which is more obvious for RMSSD than for SDNN. The RMSE in all 40 sets of data is 4.91 ms for SDNN and 9.10 ms for RMSSD. The HRV estimation results of 10 experimental subjects are listed in Table 1. One set of sampling data was selected from four sets for each subject. The estimated results of MEAN obtained by the proposed algorithm are very close to the true values, with the estimation errors of all 10 subjects being less than 3 ms. For SDNN, 9 of 10 experimental subjects' estimation errors are within 5 ms, and the maximum error is 5.4 ms. For RMSSD, the estimation errors of 9 subjects are within 10 ms, and the maximum estimation error is 10.9 ms.

Discussion
As an important method for the noncontact detection of human physiological parameters, biological radar collects the vibrations outside the human heart from RF signals, analyzes the phase changes, extracts heartbeats, and estimates the corresponding parameters through signal processing algorithms. In this article, a set of algorithms is proposed for HRV monitoring via FMCW radar. The hardware platform we use actually has a sampling frequency of 100 Hz, which means that the interval between two adjacent sampling points in the heartbeat waveform is 10 ms. The experimental results show that the measurement errors of the proposed algorithm for the HRV indexes are less than 10 ms, so we believe that the accuracy is satisfactory.
The measurement of the time intervals of heartbeats based on RF methods depends largely on the quality of echo signals and a high signal-to-noise ratio (SNR) of the heartbeat signal. Radar equipment with better performance and measurement constraints, such as the position relative to the radar and the body postures of the subject, are usually required to obtain higher-quality signals. Locating the radar close to the heart usually leads to higher signal quality. In some systems, such as the PhysioChair [40], the stress assessment method [41], and the method for monitoring a driver's heartbeat information [42], the radar is installed in the position of the seat back directly opposite to the chest, and obvious heartbeat signals can be seen in the chest movement waveform demodulated from radar echoes. A simple processing method like a bandpass filter is sufficient to extract heartbeat signals with high quality. However, it is difficult to objectively evaluate the effects of various HRV monitoring methods based on millimeter-wave radar with different RF equipment and in different monitoring conditions. Therefore, we selected several relevant studies that used the same or very similar RF devices as we used and compare these studies in Table 2, where, 'Distance' represents the distance between the antenna and the subject, and 'HR ACC' represents the accuracy of heart rate measurement. Compared with the existing common HRV extraction algorithms, such as the heartbeat extraction algorithm based on decomposition [46], which has an RMSE of MEAN, SDNN, and RMSSD of 4.65 ms, 10.31 ms, and 9.18 ms, respectively, the proposed algorithm achieves higher accuracy, with the RMSE being 3.24 ms, 4.91 ms, and 9.10 ms. Although the stochastic resonance algorithm [30] can estimate SDNN with a relative error of less than 6.53%, it is not ideal for some HRV indices (such as pNN50), as it cannot accurately measure the time duration of each heartbeat. The algorithm based on Gaussian pulse train modeling and frequency-time phase regression (FTPR) [34] has a similar problem.
A bandpass filter (BPF) is a very common method to extract heartbeat signals from radar signals. It is also used in HRV monitoring based on millimeter-wave radar [23,34]. However, the extraction effect is significantly affected by the quality of the original data. To compare the proposed algorithm with the algorithm based on the bandpass filter, a 0.8-1.5 Hz IIR bandpass filter was designed to process the same 40 sets of data. The measurement errors between MEAN, SDNN, and RMSSD obtained from each group of data and the gold standard were calculated. The CDF of errors is shown in Figure 11a-c. It can be seen that the accuracy of the HRV index calculated by the proposed algorithm is significantly higher than that calculated by IBIs extracted by the bandpass filter. For the MEAN index, the measurement error of 90% of the experimental data is less than 5 ms, and the maximum error is less than 15 ms. However, only 15% of the experimental data with bandpass filters have measurement errors of less than 5 ms, and even more than 20% of the experimental data have measurement error greater than 50 ms. For the SDNN index and the RMSSD index, the proposed algorithm can also achieve measurement errors of less than 13 ms and 24 ms, but for the bandpass filter, nearly half (47.5%) of the experimental data have SDNN index measurement errors exceeding 50 ms. Only 37.5% of the experimental data have RMSSD index measurement errors within 50 ms. It can be concluded that the HRV index calculated based on the heartbeat signal extracted by a simple bandpass filter has a large error, and its accuracy cannot meet the monitoring needs of the HRV index considering that the normal values of SDNN and RMSSD of normal people are usually tens of milliseconds. Calculating the acceleration of the chest wall motion signal to enhance the heartbeat signal is a key part of the proposed algorithm. A similar method was also adopted in [38], showing a good ability to suppress the harmonics of breath and the potential for HRV monitoring. However, the HRV of only one set of data was analyzed in their experiment, and the relative error was more than 30% for RMSSD. It should be pointed out that IBIs with high accuracy cannot be extracted only by the smoothed acceleration signal. That is, the optimization of segmentation to adjust the time interval of each heartbeat in the proposed algorithm is also necessary. To illustrate this point, the HRV index calculated by IBIs extracted only by differential strengthening and smoothing is also shown in Figure 11a-c. It can be seen that for the MEAN index, the difference between differential enhancement and the proposed algorithm is not obvious, because the position adjustment of segmentation points between each heartbeat does not affect the total duration of all heartbeats. Both methods can achieve measurement errors of less than 15 ms. However, for SDNN and RMSSD, which are more sensitive to the measurement error of each IBI, the proposed algorithm is significantly better than that without the optimization of segmentation. For SDNN, the maximum measurement error obtained by only differential enhancement exceeds 45 ms, while for RMSSD, more than 10% of the results obtained by only differential enhancement have measurement errors exceeding 50 ms. The use of only differential enhancement cannot meet the monitoring accuracy requirements of HRV indexes either, especially for SDNN and RMSSD indexes. Table 3 quantitatively shows the comparison between the proposed algorithm and the above two algorithms in terms of accuracy, where BPF represents the algorithm based on the bandpass filter, and DE represents the algorithm with only differential strengthening without the optimization of segmentation. The algorithm in [39] achieved a good measurement accuracy for IBI measurement with an average error of 3.2 ms. However, the amount of computation required by this algorithm is very large, making it difficult to monitor in real time. It took several minutes to run the corresponding script on an ordinary PC to process a set of 15 s sampling data. The computational complexity problem also exists in the algorithm in [35]. It achieved an average error of 12.2 ms for IBI measurement and successfully realized myocardial infarction detection based on HRV. However, this algorithm contains a lot of mode decomposition and optimal segmentation. The proposed algorithm achieves a comparable measurement error with much less computation. The average measurement error of the proposed algorithm is 9.3 ms, and it only takes a few seconds to complete a set of 15 s sampling data by running the script of the proposed algorithm under the same conditions. The proposed algorithm achieves a better balance between accuracy and computational complexity, which shows the potential of real-time monitoring. This study also reveals the potential of the proposed algorithm for cardiovascular system diagnosis. We noticed that most of these products claim that they can measure the heart rate with less than 2% error. Table 2 shows that the accuracy of our algorithm for heart rate measurement is higher than 98%, indicating that the proposed algorithm has the potential for clinical application in terms of accuracy. It should be pointed out that the above accuracy was obtained on the basis of data measured in a laboratory environment, which is relatively simple, and the number of experimental subjects was relatively small. A larger range of experimental verification is required before clinical application.
The proposed method also has some limitations. Firstly, its anti-interference capability is still limited in some cases, such as when monitoring multiple targets at the same time. Figure 12 shows the IBIs of one target when there are two individuals in the field of view (FOV) of the radar. It can be seen that the proposed algorithm does not extract the IBIs very well in this case. This may be due to the superimposition of chest wall motion signals from two individuals. Although MVDR beam-forming technology was adopted in our algorithm to suppress interference from other directions, because the number of antennas in our hardware equipment is limited, the echo signal from another individual is not completely suppressed when extracting the heartbeat information of one of the monitored targets. In addition, it was found in the experiments that when subjects cough or scratch, large interference will be introduced into the acceleration signal, leading to a significant deviation of the IBI estimation results from the real value. We also noticed that for a long period of data (e.g., 5 min), the measurement effect of this algorithm for IBIs is not ideal, because the acceleration waveform of a single heartbeat will also change over time. Therefore, it is difficult to find a unified template signal to segment the long period of acceleration.

Conclusions
This paper proposes an algorithm for HRV monitoring. Firstly, it extracts the echo signal of the heart through the distance and azimuth angle relative to the antenna array parameters, which are measured before the experiments. The acceleration of the chest wall motion signal is calculated to suppress respiration and strengthen the heartbeat. Then, the acceleration is smoothed by short-time average power, from which the cut-off points are estimated. Finally, IBIs are extracted by a joint optimization algorithm, and HRV features are estimated. Experimental results show that this algorithm can achieve good measurement accuracy in a laboratory environment. However the measurement performance degrades in a complex real scenario or a long-time monitoring scenario. In an environment that is not suitable for ECG and PPG measurement, the measurement method based on RF technology can provide a new approach for the measurement and analysis of HRV. In the future, we plan to apply the proposed algorithm to sleep monitoring. The trunk of the subject is relatively still during sleep time, so the HRV indicators can be estimated more accurately. We wish to help detect some diseases, such as OSA, by monitoring the HRV indicators of subjects during sleep.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Data Availability Statement: Data will be made available on request.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: